/* Copyright (c) 2010, Carl Burch. License information is located in the
 * com.cburch.logisim.Main source code and at www.cburch.com/logisim/. */

package com.cburch.draw.shapes;

import com.cburch.logisim.data.Bounds;

public class CurveUtil {
    private CurveUtil() {
    	
    }

    /**
     * getBounds and findNearestPoint are based translated from the ActionScript
     * of Olivier Besson's Bezier class for collision detection. Code from:
     *   http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
     */

    // a value we consider "small enough" to equal it to zero:
    // (this is used for double solutions in 2nd or 3d degree equation)
    private static final double zeroMax = 0.0000001;

    // note: p0 and p2 are endpoints, p1 is control point
    public static Bounds getBounds(double[] p0, double[] p1, double[] p2) {
        double[] A = computeA(p0, p1);
        double[] B = computeB(p0, p1, p2);

        // rough evaluation of bounds:
        double xMin = Math.min(p0[0], Math.min(p1[0], p2[0]));
        double xMax = Math.max(p0[0], Math.max(p1[0], p2[0]));
        double yMin = Math.min(p0[1], Math.min(p1[1], p2[1]));
        double yMax = Math.max(p0[1], Math.max(p1[1], p2[1]));

        // more accurate evaluation:
        // see Andree Michelle for a faster but less readable method
        if (xMin == p1[0] || xMax == p1[0]) {
            // u where getTan(u)[0] == 0
            double u = -A[0] / B[0];
            u = (1 - u) * (1 - u) * p0[0] + 2 * u * (1 - u) * p1[0] + u * u * p2[0];
            if (xMin == p1[0]) {
                xMin = u;
            } else {
                xMax = u;
            }

        }
        if (yMin == p1[1] || yMax == p1[1]) {
            // u where getTan(u)[1] == 0
            double u = -A[1] / B[1];
            u = (1 - u) * (1 - u) * p0[1] + 2 * u * (1 - u) * p1[1] + u * u * p2[1];
            if (yMin == p1[1]) {
                yMin = u;
            } else {
                yMax = u;
            }

        }

        int x = (int) xMin;
        int y = (int) yMin;
        int w = (int) Math.ceil(xMax) - x;
        int h = (int) Math.ceil(yMax) - y;
        return Bounds.create(x, y, w, h);
    }

    private static double[] computeA(double[] p0, double[] p1) {
        return new double[] { p1[0] - p0[0], p1[1] - p0[1] };
    }

    private static double[] computeB(double[] p0, double[] p1, double[] p2) {
        return new double[] {
                p0[0] - 2 * p1[0] + p2[0],
                p0[1] - 2 * p1[1] + p2[1] };
    }

    // returns t:Number, pos:Point, dist:Number, nor:Point 
    // (costs about 80 multiplications+additions)
    // note: p0 and p2 are endpoints, p1 is control point
    public static double[] findNearestPoint(double[] q,
            double[] p0, double[] p1, double[] p2) {
        double[] A = computeA(p0, p1);
        double[] B = computeB(p0, p1, p2);

        // a temporary util vect = p0 - (x,y)
        double[] pos = { p0[0] - q[0], p0[1] - q[1] };
        // search points P of bezier curve with PM.(dP / dt) = 0
        // a calculus leads to a 3d degree equation :
        double a = B[0] * B[0] + B[1] * B[1];
        double b = 3 * (A[0] * B[0] + A[1] * B[1]);
        double c = 2 * (A[0] * A[0] + A[1] * A[1]) + pos[0] * B[0] + pos[1] * B[1];
        double d = pos[0] * A[0] + pos[1] * A[1];
        double[] roots = solveCubic(a, b, c, d);
        if (roots == null) {
            return null;
        }


        // find the closest point:
        double tMin = Double.MAX_VALUE;
        double dist2Min = Double.MAX_VALUE;
        double[] posMin = new double[2];
        for (double root : roots) {
            double t;
            if (root < 0) {
                t = 0;
            } else if (root <= 1) {
                t = root;
            } else {
                t = 1;
            }

            getPos(pos, t, p0, p1, p2);
            double lx = q[0] - pos[0];
            double ly = q[1] - pos[1];
            double dist2 = lx * lx + ly * ly;
            if (dist2 < dist2Min) {
                // minimum found!
                tMin = root;
                dist2Min = dist2;
                posMin[0] = pos[0];
                posMin[1] = pos[1];
            }
        }

        if (tMin == Double.MAX_VALUE) {
            return null;
        } else {
            return posMin;
        }
    }

    private static void getPos(double[] result, double t,
            double[] p0, double[] p1, double[] p2) {
        double a = (1 - t) * (1 - t);
        double b = 2 * t * (1 - t);
        double c = t * t;
        result[0] = a * p0[0] + b * p1[0] + c * p2[0];
        result[1] = a * p0[1] + b * p1[1] + c * p2[1];
    }

    // a local duplicate & optimized version of com.gludion.utils.MathUtils.thirdDegreeEquation(a,b,c,d):Object
    //WARNING: s2, s3 may be non - null if count = 1.
    // use only result["s"+i] where i <= count
    private static double[] solveCubic(double a, double b, double c, double d) {
        if (Math.abs(a) > zeroMax) {
            // let's adopt form: x3 + ax2 + bx + d = 0
            // multi-purpose util variable
            double z = a;
            a = b / z;
            b = c / z;
            c = d / z;
            // we solve using Cardan formula: http://fr.wikipedia.org/wiki/M%C3%A9thode_de_Cardan
            double p = b - a * a / 3;
            double q = a * (2 * a * a - 9 * b) / 27 + c;
            double p3 = p * p * p;
            double D = q * q + 4 * p3 / 27;
            double offset = -a / 3;
            if (D > zeroMax) {
                // D positive
                z = Math.sqrt(D);
                double u = ( -q + z) / 2;
                double v = ( -q - z) / 2;
                u = (u >= 0)? Math.pow(u, 1. / 3) : -Math.pow( -u, 1. / 3);
                v = (v >= 0)? Math.pow(v, 1. / 3) : -Math.pow( -v, 1. / 3);
                return new double[] { u + v + offset };
            } else if (D < -zeroMax) {
                // D negative
                double u = 2 * Math.sqrt( -p / 3);
                double v = Math.acos( -Math.sqrt( -27 / p3) * q / 2) / 3;
                return new double[] {
                    u * Math.cos(v) + offset,
                    u * Math.cos(v + 2 * Math.PI / 3) + offset,
                    u * Math.cos(v + 4 * Math.PI / 3) + offset };
            } else {
                // D zero
                double u;
                if (q < 0) {
                    u = Math.pow( -q / 2, 1. / 3);
                } else {
                    u = -Math.pow( q / 2, 1. / 3);
                }

                return new double[] {
                    2*u + offset,
                    -u + offset };
            }
        } else if (Math.abs(b) > zeroMax) {
            // a equals 0 then actually a 2nd degree equation:
            // form : ax2 + bx + c = 0;
            a = b;
            b = c;
            c = d;
            double D = b*b - 4*a*c;
            if (D <= -zeroMax) {
                // D negative
                return null;
            } else if (D > zeroMax) {
                // D positive
                D = Math.sqrt(D);
                return new double[] {
                        ( -b - D) / (2 * a),
                        ( -b + D) / (2 * a) };
            } else {
                // D zero
                return new double[] { -b / (2 * a) };
            }
        } else if (Math.abs(c) > zeroMax) {
            // a and b are both 0 - we're looking at a linear equation
            return new double[] { -d / c };
        } else {
            // a, b, and c are all 0 - this is a constant equation
            return null;
        }
    }

    // Translated from ActionScript written by Jim Armstrong, at
    // www.algorithmist.net. ActionScript is (c) 2006-2007, Jim Armstrong.
    // All rights reserved.
    //
    // This software program is supplied 'as is' without any warranty, express,
    // implied, or otherwise, including without limitation all warranties of
    // merchantability or fitness for a particular purpose.  Jim Armstrong shall
    // not be liable for any special incidental, or consequential damages,
    // including, without limitation, lost revenues, lost profits, or loss of
    // prospective economic advantage, resulting from the use or misuse of this
    // software program.
    public static double[] interpolate(double[] end0, double[] end1, double[] mid) {
        double dx = mid[0] - end0[0];
        double dy = mid[1] - end0[1];
        double d0 = Math.sqrt(dx * dx + dy * dy);

        dx = mid[0] - end1[0];
        dy = mid[1] - end1[1];
        double d1 = Math.sqrt(dx * dx + dy * dy);

        if (d0 < zeroMax || d1 < zeroMax) {
            return new double[] { (end0[0] + end1[0]) / 2,
                    (end0[1] + end1[1]) / 2 };
        }

        double t = d0 / (d0 + d1);
        double u = 1.0 - t;
        double t2 = t * t;
        double u2 = u * u;
        double den = 2 * t * u;

        double xNum = mid[0] - u2 * end0[0] - t2 * end1[0];
        double yNum = mid[1] - u2 * end0[1] - t2 * end1[1];
        return new double[] { xNum / den, yNum / den };
    }
}